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Abstract 



Linear differential equations of arbitrary order with polynomial coefficients are 
i-^ . considered. Specifically, necessary and sufficient conditions for the existence of poly- 

nomial solutions of a given degree are obtained for these equations. An algorithm 
' to determine these conditions and to construct the polynomial solutions is given. 

■ The effectiveness of this algorithmic approach is illustrated by applying it to several 

differential equations that arise in mathematical physics. 
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1 Introduction 



Differential equations of the form Y2k=o ^kV^'^'' = where is a polynomial of degree < k 
{I < k < N) have been studied by many authors, notably Bochner [1] and Brenke [3] 
^ ! for N = 2, and Krall [9] and Littlejohn [11] (see also [12]) for orthogonal polynomial 

■ solutions. However, the case when the polynomials ak have arbitrary degree has not been 

investigated as extensively. In their recent paper [6j , Ciftci et al. considered certain types 
of such equations that arise in mathematical physics. They specifically gave conditions 
for the existence of polynomial solutions using, in particular, the asymptotic iteration 
method (AIM) they introduced in their earlier work [5]. 

In this paper we consider linear differential equations of arbitrary order with polynomial 

coefficients of arbitrary degree. Our approach is based on linear algebra and provides not 

only a necessary and sufficient condition for the existence of polynomial solutions of such 

equations, but also an algorithmic procedure for the verification of this condition as well 

as for constructing these solutions. This is discussed in detail in Section 2. In Section 3, 

we include a Maple program that can be implemented to determine the conditions that 

guarantee the existence of polynomial solutions of any degree and to find them, depending 
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of course on the available computational power. We hope that the computer program given 
here will be useful for researchers in solving linear differential equations with polynomial 
coefficients. To illustrate the efficiency of our algorithmic procedure, the Maple code 
is implemented in several examples, namely the one- dimensional Schrodinger equation, 
planar Coulomb diamagnetic problem, Bohr Hamiltonian with Davidson potential, and 
radial Schrodinger equation with shifted potential. We point out that in |6| the authors 
stated that finding polynomial solutions is a problem that needs to be investigated and 
suggested AIM for that. 

2 Polynomial Solutions 

Throughout, P is the space of all real polynomials and P„ is the subspace of polynomials 



where D is the usual differential operator and Pfc(x) = 'YlihxdPkh^^ is a polynomial of de- 
gree dk ( with the convention that the zero polynomial has degree — oo and that D^y = y). 
Our objective is to find a necessary and sufficient condition for the equation Ly = to 
have non-trivial polynomial solutions. Although this can be achieved, for each specific 
case, by comparing coefficients (see for example the determinantal necessary condition in 
the recent interesting paper [6] on Heun's equations), or by using the Asymptotic Itera- 
tion Method in the case of second-order equations [5] , we feel that a systematic approach 
that works for differential equations of all orders and that can easily be implemented in 
a computer algebra system is more desirable. 

Assume first that for some i {0 < i < N), di > i. Let m = max {di — i) and put y = D^z. 
In this way, the equation = is equivalent to Hz = where H is the linear operator 



y ak{x)D , and Ly = has a polynomial solution of degree n > if and only if Hz = 

k=l 

has a polynomial solution of degree n + m. Clearly, for each nonnegative integer n, P„ 
is if- invariant, and H has thus the advantage over L of being directly amenable to an 
eigenvalue analysis as demonstrated below. We note here that the easier case when di < i 
for all i can be discussed almost verbatim, with obvious modifications. 

Let Ofc {k > 1) be the sequence of polynomials defined by = if A; < m and = pk-m 
if /c > m. Put afc(x) = J2h>o'^khX^, where akh = if k < h. Since, for each nonnegative 
integer n, is a scalar multiple of x" plus lower order terms, we see that the matrix 

representation of H, with respect to the standard basis Bn = {l,x, of P„ is upper 




m+N 
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triangular and its eigenvalues are the coefficients of in H{x"'). More specifically, the 
(n + l) X {n+ 1) matrix An of H operating on P„ has {i, j)-th entry ak,k+i-j{3 ~ 

k>l 

i.e. 



A 



n — ^2 ^k,k+i-jU — k)k 

-''^^ -I l<i,j<n+l 

where {j — k)k = {j — — '2) ■ ■ ■ {j — k), and where each row and column has at most 
{N + m + 1) nonzero entries. Clearly, the first m columns of An are zero and An+i is 
obtained by An by adding one row and one column at the end. As diagonal entries of An, 

n 

all the eigenvalues of the operator H are real and are given by A„ = n! > -. —7 for 

n > 1 (note that Aq = Ai = • • • = Xm-i — 0). Each eigenvalue A„ has an eigenpolynomial 
Unix) — Uno + Vnix + • • • + i/nnx"' of degree at most n and whose vector representation 
iUno, ■■■,ynn)^ in the standard basis Bn can be directly computed from the homogeneous 
upper triangular system {An — A„/)(y„o, ■■■,ynn)^ — 0. Our problem is to find necessary 
and sufficient conditions for which the operator H has an eigenpolynomial of degree n + m 
corresponding to Xn+m — 0, that is necessary and sufficient conditions for the homoge- 
neous system A„+^(y„+^,o, ■■■:yn+m,n+mV = tO have a solution (yn+m,0, ■-, Z/n+m.n+m)^ 

with yn+m,n+m — 1- This will follow from 

Lemma 1. Let A be an m x n matrix. Then the homogeneous system AX — has a 
solution X — {xi,X2, Xn)'^ with Xk ^ for some k if and only if rank{A) — rank{Ak) 
where {A^) is the matrix obtained from A by deleting the A;*'* column. 

Proof. Put A — [cjj]i<ij<„ and let be the k*^ column of A. Clearly, A and the 
augmented matrix [AfeXfe] have the same rank. Hence, 

rank{A) — rank{Ak) <^ rank[Ak:ck\ — rank{Ak) 44> the systemA^X = —Ck is consistent 
4^ there exists a solution X — {xi, • • • , Xk-i, 1, Xk+i, • • • , Xn)'^ to the system AX = 0. □ 

Since An+m-i is obtained from An+m by deleting the last column, the above lemma im- 
mediately yields that the differential equation Ly = has a polynomial solution of degree 
n > if and only if rank{An+m) — i^^nklAn+m-i)- In this case, since An+m is upper tri- 
angular, the last entry of An+m is zero i.e. Xn+m — o-kkiTn + n — k)k = 0, and therefore 

fe>i 

the last row of An+m is zero. 
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Now let Mn and be, respectively, the matrices obtained from An+m and An+m-i by 
deleting the first m zero columns. Clearly rank{An+m) = rank{An+m-i) if and only if 
rank{Mn) = rank{M'^. It is easy to see that the entry of the (n + m+ 1) x (n + 1) 

matrix M„ is ^ at+m,t+i~j{j -t)t+m = ^Pt,t+i-j{j -t)t+m- This proves the main result 

t=0 t=0 

of this note: 

j-t 

Proposition 2. Let M„ be the (n+m+1) x (n+1) matrix with {i,jy'^ entry pt, t+i- j (j ~ 

t=o 

t)t+m and let be the matrix obtained from M„ by deleting the last column. Then the 
differential equation Ly = has a polynomial solution of degree n > if and only if 
rank{Mn) = rank{M^).a 

It thus follows that if the equation Ly = has a polynomial solution of degree n > 0, then 
Pt,t+m{n — i)(t+m) — 0, and since has n columns, rank{Mn) — rank{M^) 



implies that rank{Mn) < n and so every (n + 1) x [n + 1) submatrix of M„ has zero 
determinant. This generalizes Theorems 5 and 6 of [B]. 

3 Maple Code and Examples 

In this section we illustrate the effectiveness of the algorithmic approach of Section 2 by 
applying it to four differential equations that appear in [6]. These arise in mathematical 
physics and, more precisely, in the study of solutions to Schrodinger equation [lOj, planar 
Coulomb diamagnetic problem Bohr Hamiltonian with Davidson potential [2] and 
radial Schrodinger equation with shifted potential [SI [71 Ej. The examples show how to 
implement the method algorithmically to determine the conditions for the existence of 
polynomial solutions and also to calculate the corresponding polynomial solutions. 

Example 1 

As a first example, we consider the linear second order ODE arising in the study of one 
dimensional Schrodinger problems . The investigation of Krylov and Robnik [10] about 
polynomial solutions of one dimensional Schrodinger problems leads to investigation of 
polynomial solutions of the following differential equation. The conditions for existence 
of polynomial solutions of this ODE have also been discussed by Ciftci et al. in the recent 




t>i 
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paper [6], by a different approach. 



X 



dx' 



y{x) + a [x^ - l) —y (x) + x + g)y {x) = 



(3.1) 



Here, we apply our method to determine existence conditions as well as to compute the 
corresponding polynomial solutions of the above differential equation. A Maple code is 
also provided to show how the method can be implemented algorithmically using software. 

For (7 = 0, the algorithmic procedure of Section 2 can be implemented to generate a 
sequence of even degree polynomial solutions of ODE (13. ip . A further analysis of these 
solutions yields the following result. 

• For (7 = 0, ODE (13. ip admits polynomial solutions of degree n = 2m (m > 1), with 

/3 = —{an + n"^ — n), given by 



-mi) 



^ij.2m—2i 



^ {a + 2n- 3)(a + 2n - 5) ■ • ■ (a + 2n - 3 - 2(z - 1)) 

It should be noted that, for (7 = 0, no lower odd degree polynomial solutions of 
ODE ([31]) were found. 

For (7 7^ some examples of polynomial solutions, found using the construction pro- 
cedure of Section 2, are given in Table 1 below. 



n 



/3 
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Polynomial solution of ODE (13. ip of degree n 



—a 



x±l 



-2a - 2 



2 ^ V2a(2a+3) 
a+2 



+ — 

~ a+2 



-3a — 6 



±Vl5a + + aA 



. 3 , 12a(9+2a+A) 

aVl5+5a+A{-3-2a+A) 
+ Ji^. . X ± 24a2 



X 



-3-2a+A' 



a(-3-2a+A)V15+5a+A 



±-\/15cfc + — aA 



Ax^ ± 



12a(-9-2a+A) 



a^/15+ba-A(■i+2a+A) 

24q -j- 24a^ 

3+2q:+A"'' "T 



X 



a(3+2a+A)V15+5a-A 

where 

A = V153 + 96a + 16a^ 



Table 1: 



In general, for any given n and a, the algorithmic procedure of Section 2 can easily be 
implemented to determine /3, g for which ODE (13. ip admits polynomial solutions of degree 

n as well as to compute the corresponding polynomial solution. We provide below a set of 
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Maple commands that can be used as template to compute polynomial solutions of given 
degree n of ODE (13 .ip for any given value of a. For illustration we take a = and look 
for a solution of degree n = 6. The following Maple code determines (3 = 15, g = 3(750)^/^ 
and computes the corresponding polynomial solution of degree 6 of ODE (13. ip as 

7 oM , VSO (-65 + 12 730) , v^TOSO (6 V30 - 30) , 

y(x) = 7x^ + — lOSO^/^x^ + 126 — ^-x^ + 420 ^-x^ 

30 66^30 - 360 66^30 - 360 

(-72 + 6^30) „ (-72 + 6^30) 750=^/^ , , 

+ 1575^ = ^-x^ + 630- — ^ =-x 3.2 

66 V30 - 360 (66 ^30 - 360) (-60 + 5 VSO) 

(-72 + 6^30) 730 

(66 V30 - 360) (-60 + 5 ^30) 

The Maple code with brief explanations is presented below. 
restart: 

with(Linear Algebra) : 
a := -15/2 : 
N ■=2: 

pcoeff := Array (0 .. N): 
pcoejf[Oj := l3x + g : 
PCoeff[l] := a(x2 - 1) : 
pcoeff[2] .-= x^ : 

These commands define the order of the ODE, the value of a and the coefficients Pk{x) 
of the operator Ylk=o ^^y- Next we determine the value of m and the coefficients a^i^x) 
of the operator Ylk^^i^ using the following commands. 

Vm := Array{Q..N) : 

dk := Array{O..N) : 

for i from to N do 

dk[i] := degree (pcoeff [i], x) 
Vm[i]:= degree(pcoeff[i] ,x)-i end do: 

m := maxiVm) : 

dkmax := max{dk) : 

acoeff := Array (1 .. m+N): 

for i from m to (m+N) do 
acoeff[i]:= pcoeff[i-m] end do: 
Matrices An+m, ^n+m-i of the procedure of Section 2 are computed by the following set 
of commands in which " soldegree" , An and Anprime respectively denote the degree of 
the sought polynomial solution, the matrix An+rn and the matrix An+m~i- 

soldegree := 6: 
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n :— soldegree+m: 

cAM :— max(dkmax, n+m+N): 

rAM m+N: 

AM Array (1 .. rAM, .. cAM): 
for i from m to rAM do 

AM[i,0]:— coeff(x acoeff[i],x^) end do: 
for i from m to rAM do 

for j from 1 to cAM do 

AM[i, jj coeff(acoeff[iJ, x^ ) 

end do 
end do 

An Matrix(n+l,n+l): 
for i from 1 to n+1 do 
for j from 1 to n+1 do 
for k from 1 to rAM do 

if (k+i-j) > and (k-j+1) < 

then AnfiJJ:^ An[i,j]+AM[k,k+H] (jS)! 
end if 
end do 
end do 
end do 

Anprime := Matrix(n+l,n) 
for i from 1 to n+1 do 

for j from 1 to n do 

Anprime[i, j] :— An[i, j] 

end do 
end do 

At this stage we have Rank{An) ^ Rank{Anprime) . Next we determine the values of 
parameters so that Rank{An) equals Rank{Anprime) . The first condition employed is 
the vanishing of the last diagonal entry of the upper triangular matrix An which deter- 
mines f3 by the following commands. ansbeta :— solve(An[n+l, n+1] — 0, j3): 
/3 := ansbeta 

Implementing the fact that the 7x7 submatrix, obtained by deleting the last zero row of 
An, must have zero determinant provides the value of g via the commands below. 

Andet := Matrix(n,n): 

for i from 1 to n do 
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for j from 1 to n do 

Andet[i, j] := An[i, j+1] 
end do 
end do 

Determinant ( Andet ) : 

ansdet := solve(Determinant(Andet) = 0, g) 
This leads to seven roots. At this stage a root needs to be chosen before checking the 
rank condition. For illustration we choose g = 3(750)^/^. 

g:= 3(750)1/^ : 

Rank(An) 

Rank( Anprime ) 

As the rank condition is satisfied so the desired polynomial solution can be obtained by 
the following set of commands. 

kern := NullSpace(An) 
The output of the above command contains two vectors kerni and kern2 with kern2 = 
(1, 0, 0, 0, 0, 0, 0)"^ so the vector kerni is used as follows to find the solution. 

Vkern := kerni : 

solm := 0: 

for i from 1 to (n+1) do 

solm:= solm + Vkern[i] x^^^ end do: 
soln := diff(solm, x) 
The output soln provides the solution given in Equation (13. 2p . 

Example 2 

Consider the ODE 

—y (x) + (p - 2a;2) —y (x) + (5 x + a) y (x) = (3.3) 

The question of investigating the polynomial solutions of ODE (13. 3p arises from the study 
of polynomial solutions of Coulomb diamagnetic problem by Chhajlany and Malnev [1]. 
Ciftci et al. |6l Eq.18,19] provide conditions for the existence of polynomial solutions of 
ODE (13. 3p . Here we implement our procedure to demonstrate how to generate polynomial 
solutions of ODE (13. 3p in a straightforward algorithmic manner. 

For general a ^ some examples of polynomial solutions listed in Table 2 below are 
obtained by adapting the Maple code presented in Example 1. The conditions on the 
parameters 6 and p, for having these solutions, are also determined. 
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p 



Polynomial solution of ODE (13. 3p of degree n 



2x — a 



0^+16 

aa 



„2 _ a _|_ 
± "2 ma— 



5af 
18 



± 



2Va"t-54a 
9 



4x^ — 2ax'^ 



^ (±5 0^+4 ^a{o3_54)a=p216jax 
±«2+2 ^a(o3-54) 



1 ±41 a5+40a3^Q {«3-54)=Fl728 a2_432 (a3_54) 
_54 =fca2+2^a (a3-54) 



J_ 5a3+192ifc3A 
64 Q 



c;™4 5a 3 i lo 



15 



hFa^±768 a3+Q,3^zp4096_64A)x2 



2 '32 a(±3o3=pi92+5A) 
5 (=Fo''Tl28o«+a''A±8192a3_i728a3^_|_262144+4096 A)x 

64 (±3a3=pi92+5^)(=Fa3±64+yl) 

-J- 5 B 

+ 2048 a2(_|_3ct3zpi92+5A)(=FQ3-|-64+^) 

where 



A = Va^ - 384 a3 + 4096 and 

B = ^ 2816 + a^A ± 524288 - 5184 

±19922944 + 200704 t 150994944 - 2359296 A. 



Table 2: 

Depending on the available computational power, for a given value of a, the algorith- 
mic procedure of Section 2 can be implemented, as in Example 1, to compute polynomial 
solutions of ODE (13. 3 p of any given degree n. As example the following solutions of degree 
9 and 25 are found. 



a = 0, n = 9 implies 6 = 18, p = ^(144830)^/'^ and the polynomial solution of 
degree 9 of ODE (13. 3p given by 



y[x) 



lOx^ - ^ ^/l44830x^ - 120x^ + - 144830^/3x5 + — 4^UmOx^ 



10314 



8478 



35 
4239 



35 

19803534 



7 



X 



144830^/^x2 + v^l44830x + 

8575 245 8575 



For a = 0, n = 25 with 5 = 50 and p = the polynomial solution of degree 25 of 
ODE (I33D is given by 



y{x) = 26x^5 - 2600x^2 + 100100 x^^ - 1901900 x^^ + 19019000 x 
- 98898800 x^° + 247247000 x^ - 247247000 x^ + 61811750 x 
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Example 3 

The analysis of solutions of the Bohr Hamiltonian for Davidson potential leads to inves- 
tigation of exact solutions of a differential equation [2], Eq.49] which can be rewritten as 
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la Eq.22] 

cf d 
x—y{x) - (2x2-2/2-2) —y{x) - {2 fi + 3 - e) xy {x) = (3.4) 
dx^ ax 

Adapting the Maple code of Example 1 for ODE (13 .4^ readily generates polynomial solu- 
tions of a given degree n. Some examples for solutions of even as well as odd degrees are 
provided below in Tables 3 and 4 respectively. 



n 


e 


Polvnomial solution of ODE (13.41) of degree n = 2m 





2/i + 3 


1 


2 


2/i + 7 


2x2 - (2/X + 3) 


4 


2/i + 11 


4x^-4 (2/i + 5)x2 + (2/i + 3)(2/i + 5) 


6 


2/i + 15 


8 x^ - 12 (2 /i + 7) x^ 

+6 (2 /i + 5) (2 + 7) x^ - (2 /i + 3) (2 /i + 5) (2 /i + 7) 


8 


2/i + 19 


16 x« - 32 (2 /i + 9) x** + 24 (2 /i + 7) (2 /i + 9) x^ 

-8 (2 /i + 5) (2 /i + 7) (2 /i + 9) x2 + (2 /i + 3) (2 + 5) (2 /i + 7) (2 /X + 9) 


10 


2/i + 23 


32 x^'J - 80 (2 /i + 11) x« + 80 (2 /i + 9) (2 /i + 11) x*^ 
-40 (2/i + 7) (2/i + 9) (2/x + ll)x^ 
+10 (2 /i + 5) (2 /i + 7) (2 /X + 9) (2 /i + 11) x^ 
-(2/i + 3)(2/x + 5)(2/i + 7)(2/i + 9) (2//+ 11) 



Table 3: 



It should be noted that for this case polynomial solutions can easily be computed with- 
out much computational cost. For instance, Maple could compute polynomial solution of 
degree n = 100 in computational time of 1.8 seconds. 



n 


e 


/i 


Polvnomial solution of ODE (13. 4p of degree n = 2m + 1 


1 


2/2 + 5 


-1 


X 


3 


2/2 + 9 


-1 


— q 2 


-2 


X'^ 


5 


2/2+13 


-1 


^5 _ _^ 


-2 




-3 


X^ 


7 


2/2 + 17 


-1 


™7 _ 21 5 I 1U5 3 _ 1U5 


-2 


x''' - 7x^ + ^x^ 


-3 


7 7 "i 


-4 


— 7 ^ 

X 



Table 4: 



A further analysis of the odd lower degree polynomial solutions of ODE ( 13. 4p . given in 
Table 4, yields the following family of solutions of degree n = 2m + 1. 



For = 0, 1, 2, 3 ODE (13.40 admits the following class of polynomial solutions of degree 
n = 2m + 1 {m> k) with /i = — (m + 1 — /c) and e = 2m + 2A; + 3. 



• If A; = (m > 0) 




• if A; = 1 (m > 1) 




2m + 1 
2 



,2m-l 



• iffc 



2,3 (m > k) 



y = x- 




k (-l)*(2m + l)(2m - 1) ■ ■ ■ (2m + 1 - 2(i - 1)) 



X' 



,2m+l-2j 



i=l 



+ 



(-l)'=(2m + l)(2m - 1) • ■ ■ (2m + 1 - 2{k - 1)) 

2^ 



Example 4 

As a final example, we consider a question related to the investigation of the radial 
Schrodinger equation with shifted Coulomb potential which has been discussed recently 
in [El [71 IE]. The anstaz of [6l Eq.35] that the radial Schrodinger equation admits a solution 
which vanishes at the origin and at infinity leads to the question of obtaining solutions of 
the following differential equation; the reader is referred to [6] for details. 



This is a particular case of the confluent Heun equation whose polynomial solutions can 
be studied algorithmically using our procedure. While discussing the question of poly- 
nomial solutions of ODE (13. 5p . Ciftci et al. [6] provide conditions on parameters a, (3 
to have polynomial solutions. In particular a table was provided which listed conditions 
on parameters for the existence of polynomial solutions for n = 1,2,3,4. However, it was 
pointed out in [6] that finding the corresponding polynomial solutions is an open problem 
that remains to be solved. For a given value K and the given degree n of the required 
polynomial solution, adapting the Maple code of example 1 for ODE (13. 5p can determine 
conditions on parameters a, /3 as well as generate corresponding polynomial solutions 
of ODE (13. 5p . In Table 5 below, we demonstrate this by providing some examples of 
polynomial solutions of ODE ( 13. Sp of degree n = 1,2, 3, 4, 5. 



x{x + (3)-^y{x) + i-2ax^ + 2{K + l- af3)x + 2f3{K + 1)) —y{x) 
dx'^ ^ ' ax 

+ {{-2a{K + 1) + 2Z)x - 2a(3{K + 1)) y{x) = 



(3.5) 
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n 


a 


K 


/3 


Polynomial solution of ODE (j3.5p of degree n 


1 


z 

K+2 


anv 


Z 




2 


Z 


any 


(3 /<+6± ^^-2+8 ii'+12j {-ft:+3) 


o 2 o iK+3){3K+6±VK-^+8K+12){2K+3)x 


iC+3 


2(X+2)Z 


(ii-+2)Z(A'±VX^+8i<-+12j 
^ (■3K+6±V-ft:2_,_8^_,.12)(_ft:+3)^(2ir+3) 

1 I 

(x+2)z2(_ft:±V-ftr24,8ft:+i2) 


3 


Z 

K+i 


3 

2 


5 

2Z 


4a;3 + fx^ 


2i> 
27, 


/|„3 1 iiU 2 I iSYb _|_ 93Y5 


4 


Z 
K+b 


3 
2 


49 
2Z 


r 4 1 1435 x'^ 1 18375 , 3109295 x , 21765065 
^•^ "T 4 Z ~'~ 2 Z2 ' 32 Z» ' 642* 


7(±15+v^) 
20Z 


Z[±1+V6E) (±1+^65)^2 


5 


Z 

K+& 


3 
2 


81 
2Z 


6x^ + 891 f + 51030 ^ + 1390932 |^ 

1 282195171 X 1 2539756539 
1 16 Z4 1 32^5 



Table 5: 
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